Fluctuations and Dispersal Rates in Population Dyanmics 
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Dispersal of species to find a more favorable habitat is important in population dynamics. Dis- 
persal rates evolve in response to the relative success of different dispersal strategies. In a simplified 
deterministic treatment (J. Dockery, V. Hutson, K. Mischaikow, et al., J. Math. Bio. 37, 61 (1998)) 
of two species which differ only in their dispersal rates the slow species always dominates. We 
demonstrate that fluctuations can change this conclusion and can lead to dominance by the fast 
species or to coexistence, depending on parameters. We discuss two different effects of fluctuations, 
and show that our results are consistent with more complex treatments that find that selected 
dispersal rates are not monotonic with the cost of migration. 

PACS numbers: 87.23.Cc, 87.10Mn 
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Dispersal plays an essential role in the dynamics of 
populations. This effect is particularly important when 
different parts of a habitat are more or less desirable 
because of a distribution of resources, varying temper- 
atures, etc. Rates of dispersal can evolve in response to 
such distributions In this paper we try to under- 
stand this evolution by considering the competition of 
two species, "fast" (F) and "slow" (S), that differ only 
by their dispersal rate. We ask who "wins" , i.e. whether 
one species drives the other to extinction in the long-time 
limit. We also consider the possibility of coexistence. 

There are two contradictory guesses that one might 
make. On the one hand, if a species is well adapted to 
part of a habitat, it would be a waste of resources to 
wander away immediately; this is exactly what species F 
does. We might postulate that F will ultimately be dis- 
placed by S provided there is some probability for S to 
find all the parts of the system. This idea is supported by 
the work of Hastings Q and Dockery, et al. Q who for- 
mulated this problem in terms of a pair of deterministic 
reaction-diffusion equations. They proved that S always 
drives F to extinction for any non-uniform habitat. 

On the other hand Hamilton and May Q considered 
a model in which F may be preferred. In this picture 
there are discrete sites that allow one organism per site. 
In each generation each organism produces offspring that 
can either remain at the site or spread at random to all 
other sites while paying a cost of migration. The adult 
in the next generation is chosen at random among the 
young present at the site. It is easy to see that S does 
not always win in this model. For example, if S has zero 
dispersal rate then S will go extinct. This results from 
the fact that there is a finite probability that a young F 
will take over any site even if it is outnumbered by S, and 
S can never return. This result depends on the stochastic 
nature of the model. Later work generalized the model 
to have more than one adult per site [5|. In this model 
there is an optimal dispersal rate which depends on the 
cost of migration. 



We will try to understand the relationship between the 
two results above by using an agent based model which is 
closely related to the reaction-diffusion equations of [3(. 
Of course, for very large populations the deterministic 
model must be correct. We explore the role of fluctua- 
tions for moderate populations. 

We are motivated to do this by the fact that stochas- 
tic effects should be important in real populations in a 
patchy environment. Further, it is known in other con- 
texts 0, 0, [H that reaction-diffusion equations can give 
misleading results. We will find that there are two ef- 
fects that lead to differences from the deterministic ap- 
proach. For a patchy environment discreteness of indi- 
viduals means that migration across unfavorable regions 
is much less frequent than that predicted by the differ- 
ential equations; this promotes coexistence. For nearly 
uniform resources there is another, more subtle, effect: 
fluctuations in population density increase death rates 
so that F, whose populations are more uniform, is fa- 
vored. To our knowledge this mechanism has not been 
discussed previously. Our simple picture explicates the 
essential features of more complex treatments based on 
data for real populations @. 

In our agent-based model the two species live at sites 
on a one-dimensional lattice with lattice constant h. Po- 
sitions are given by x j — (j — l/2)h,j — 1 . . . N x . (We can 
think of h as being comparable to the minimum, length 
scale over which the environment varies.) The numbers 
of agents at site j are F(xj,t),S(xj,t). The agents are 
identical except for their rates to jump to adjacent sites: 
Df > D s . The birth rates per agent of either type are 
non-uniform to reflect the non-uniform habitat, and are 
given by a specified function of position a(xj). To limit 
populations we assume the death rate is population de- 
pendent; we use the logistic model, so that the death 
rate per agent is (3(F + S); competition is encoded in 
the death rate. Time is in units of the inverse growth 
rate for the population in the most favorable region of 
the habitat. If F drives S to extinction, we take this as 
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an indication that in a natural population the dispersal 
rate will evolve upwards, and vice versa. 

The algorithm for our agent-based model is quite 
straightforward. First, we pick a time step dt. At each 
point there is a number of individuals of each type at 
each site. We take the probability to move left or right 
to be Dfdt for the fast species and D s dt for the slow. In 
most of the simulations we take Df = l,D s = 0.25. The 
number of agents that move in each direction is given by 
the binomial distribution. Similarly, the probability of 
birth per agent is adt. The number of births is a bino- 
mial deviate; this number is added to the original number 
at the site for each species. For deaths the corresponding 
probability for each agent is (3(F + S)dt. 

We adopt a simple landscape for a(xj) which has 
two 'oases' and a 'desert' between them: a(x) = 1 + 
7] cos(2irx / L) where L = N x h is the length of the system: 
see Figure 1(a). The parameter rj gives the range of vari- 
ation of the birth rate. If rj = the birth rate is the same 
everywhere; large rj indicates a patchy landscape and a 
desert which imposes a large cost of migration. We take 
[3 = 1/N a so that in the absence of diffusion, and for one 
species, the population at each site saturates at the car- 
rying capacity: a/fl — N (l + rj cos(2ttx / L)) . The initial 
condition is that S is localized on the left oasis, and F 
on the right, and the boundary condition is that no flux 
passes through x = or x = L. 

The result of many simulations of our model in this 
landscape is shown in Figure 1(b). The striking feature 
is that at small population densities F always wins. Fur- 
ther, for large rj and relatively small N a extinction takes 
a very long time so that in a real situation there would 
appear to be stable coexistence. Thus the region where 
F persists for very long times is non-monotonic. 

The reaction-diffusion equation version of the model 
gives different results [2, The equations are: 
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FIG. 1: (a). The landscape for a system of length 50. 
(b). Simulations results from the agent-based model; Df = 
1,D S — 0.25. Above the line marked with circles (o) S wins, 
and below, F. For large rj there is long-lived coexistence: be- 
low the line marked with squares (□) the extinction of F takes 
more than 10 times the prediction of Eq. (1) to fall below one 
individual per cell. 
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Here the parameters of the spatially discrete individual- 
based model are replaced by continuum quantities: / = 
F/h,s = S/h,D f>s = D f , s h 2 ,l3 = 0h. We should expect 
these equations to give the same results as the individual- 
based model for large N a and sufficiently smooth a(x). 

The results of Q are that for any well-behaved a(x) 
the slow species always drives the fast one to extinction, 
which is quite unlike the results of Figure 1(b). This 
result is independent of initial conditions and the detailed 
form of a. The qualitative reason for this result is that 
F cannot exploit the oases as well as S since it is always 
wandering into the deserts. However, our results show 
that even for reasonably large values of / and s effects left 
out of the partial differential equations (e.g population 



fluctuations) change the results in a qualitative way. For 
very large N a , we recover the continuum results. 

Actual populations could well be in the regime where 
the individual-based model differs from the continuum 
one, so it is useful to enquire about the sources of the 
complex behaviour. Situations where reaction-diffusion 
equations fail to capture qualitative features of stochastic 
dynamics are well known @, 0, H, HI 



11[ . For example, in 



the DLA model [8|, [12J , fluctuations in growth rate drive 
instabilities and produce fractal shapes. For Fisher-KPP 
dynamics [13|, [l4[ the discreteness of agents in an agent- 
based model dominates front propagation 

0,0, El- Dis- 
creteness also can qualitatively change low-dimensional 
dynamical models [16|. In our case both fluctuation and 
discreteness effects are important. 

First consider the case of large rj. Here the dominant 
effect is the cost of migration in crossing the desert in 
our landscape. The reaction-diffusion equations, above, 
predict an exponentially small tail of the slow species sur- 
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FIG. 2: The phase boundary marking the critical value of 
1/e (i.e., the critical value of N D for e = 1) as a function of 
r\ which separates the coexistence phase (below the line) and 
the victory of the slow species (above the line). Note the 
qualitative similarity to Figure lb, right side. 



viving into the fast oasis which eventually starts growth 
there. For discrete agents, however, this is an overesti- 
mate: we should not consider values of / which corre- 
spond to far less than one agent. Thus the tail is sup- 
pressed, and the number of agents that can cross the 
desert is very small, and are easily killed off by fluctu- 
ations. We can qualitatively represent this effect at the 
continuum level [6| by changing the birth term in the 
reaction-diffusion equations: 

a(x)f -> a(x)f@(f-e) a(x)s -> a(x)sQ(s-e), (2) 

where the function is for negative arguments and 1 
for positive, and e is a parameter of order unity. This 
means that when the density drops to very low values 
(less than one individual in a few lattice sites) , we do 
not expect to have reproduction. Numerical solutions of 
the new reaction-diffusion equations are shown in Figure 
2. Now we have a region of stable coexistence. This 
is qualitatively similar to the right side of Figure 1(b). 
In the full agent-based model, the coexistence region is 
metastable: after a long time the oasis of the slow agents 
is invaded by the fast whose transition across the desert 
is easier due to the large value of Df. 

The role of discreteness and fluctuations for small 77, 
i.e. a nearly uniform environment, is more subtle, and 
is a separate mechanism. We will argue that density 
fluctuations increase the death rate: since F smooths out 
fluctuations by diffusion, S is at a disadvantage. The 
basic reason for this is that the birth rate is linear in our 
model and the death rate non-linear. Consider the case of 
a constant, and suppose the density of the slow species 
is given by s(x) = s + Skcos(kx), where the fc's are 
chosen to satisfy no-flux boundary conditions. The Sk are 
the amplitudes of the fluctuations of various wavelengths. 
Now it is clear that the total number of births in time dt 



is independent of the fluctuations: 

dt J asdx — dtaLs . (3) 

However the total number of deaths is: 

dt J (3s 2 dx = dt(3L(s 2 + ^4/2). (4) 

Fluctuations increase the death rate. As we have noted, 
we expect the fluctuations of f(x) to be smaller, so the 
slow species is driven to extinction. 

To treat this effect in more detail we should write down 
the master equation for the stochastic process described 
above. Then we could find the differential equations for 
the averages of /, s as moments. The reaction-diffusion 
equations, Eq. (Q} should be thought of as the mean-field 
decoupling of these equations. 

For our purposes it is sufficient to average over position 
as well as over internal fluctuations. We denote these 
processes by (•). The diffusion terms average to zero, 
and the remaining terms in the moment equations must 
be of the form: 

(/) = {af) -P((f) + {fs)) 

(a) = (as)-(3((s 2 ) + (fs)). (5) 

Now we can subtract the two equations and find: 

(/-*> = (f-s}[(a) -/3{f + s}} 
+ [Cov{a, f) - Cov(a, s)] 
- /3[Var(/) - Var(s)]. (6) 

Here Cov(a, /) = (af) — (a)(f) is the covariance, and 
Var(/) = (f 2 ) — (f) 2 is the variance. 

Suppose / and s are very close to one another. Then 
the first term in the equation is zero. The second term 
favours the slow species since it can better follow the 
variation of a as in the deterministic treatment [|| . The 
third term favours the fast species since its variance is 
smaller. We can estimate the second term by setting 
either density to be of order N a a. The difference in the 
covariances is of order N {a 2 ) oc N r) 2 . In the last term 
the variances are of order N a so that /3[Var(/) — Var(s)] 
is of order unity. The two terms are equal on the division 
line on Figure 1(b), so that N a oc l/i] 2 defines the line. 
This estimate is verified by our results, see Figure [3l 

An example of these mechanisms in a more realistic 
context may be found in 9]. This is a spatially explicit 
study of the evolution of dispersal rates in a patchy en- 
vironment which is closely based on the ecology of the 
checkerspot butterfly. The model is quite complex, and 
has 15 parameters that are deduced from field observa- 
tions. In the study there is a cost of migration between 
patches of habitat, and the dispersal rate is allowed to 
evolve. As the cost of migration increases, dispersal rates 
decrease, but for large costs (analogous to large r\ in our 
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that £(100, 0.75) « 0.25 and L>(1000,0.75) « 0.15. 
The results of [3j may be interpreted to mean that 
limjVo-too D(N a , rj) = for all 77. 
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FIG. 3: A comparison of the division line between fast and 
slow dominance (the left side of Figure lb) with the estimate 
based on Eq. (|6]). The simulations are for a system of length 
10. 



treatment) there is non-monotonic behaviour, and the 
dispersal rate increases again. We propose that the un- 
derlying mechanisms in the model are revealed in our 
schematic treatment, where the evolutionary result of one 
or another species going extinct is a proxy for evolution 
of dispersal rate. We think that this kind of phenomenon, 
where fluctuation-driven effects that occur at relatively 
low densities may be operative in real systems, and will 
repay further study. In particular, an experimental real- 
ization of the model is probably possible. For example, 
In a predator-prey system of bacteria 17j fluctuation ef- 
fects on extinction have been demonstrated in dispersed 
environments. It would be very interesting to do a similar 
study in the setting of competition for resources. 

Finally we return to the question of evolution in our 
model. In 0, [B], Q there is an optimal dispersal rate. 
We have numerical evidence that this is also true for 
our model. That is, for each N a and 77 there is a 
D(N ,i]) which wins any competition and presumably 
is the evolutionarily stable state. For example, we find 
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